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Finite-Difference and Pseudospectral Time-Domain 
Methods Applied to Backwards- Wave Metamaterials 

Michael W. Feise, John B. Schneider, Member, IEEE, Peter J. Bevelacqua 



Abstract — Backwards-wave (BW) materials that have simul- 
taneously negative real parts of their electric permittivity and 
magnetic permeability can support waves where phase and 
power propagation occur in opposite directions. These materials 
were predicted to have many unusual electromagnetic properties, 
among them ampliflcation of the near-held of a point source, 
which could lead to the perfect reconstruction of the source field 
in an image [J. Pendry, Phys. Rev. Lett. 85, 3966 (2000)]. Often 
systems containing BW materials are simulated using the finite- 
difference time-domain technique. We show that this technique 
suffers from a numerical artifact due to its staggered grid that 
makes its use in simulations involving BW materials problematic. 
The pseudospectral time-domain technique, on the other hand, 
uses a collocated grid and is free of this artifact. It is also 
shown that when modeling the dispersive BW material, the 
linear frequency approximation method introduces error that 
affects the frequency of vanishing reflection, while the auxiliary 
differential equation, the Z transform, and the bilinear frequency 
approximation method produce vanishing reflection at the correct 
frequency. The case of vanishing reflection is of particular interest 
for held reconstruction in imaging applications. 

Index Terms — backwards-wave material, left-handed material, 
double-negative material, metamaterial, FDTD methods, pseu- 
dospectral time-domain method 



I. Introduction 

ELECTROMAGNETIC meta-materials with simultane- 
ously negative electric permittivity and magnetic perme- 
ability have recently received much attention. These materials 
can support a backwards-traveling wave where the phase prop- 
agation is antiparallel to the direction of energy flow. These 
materials have been identified by several names including 
left-handed or double-negative material but we shall refer to 
them here as backwards-wave (BW) materials. Their properties 
were first considered theoretically by Veselago [1] during the 
1960's but they have only been fabricated recently [2], [3]. 
They are predicted to exhibit many unusual properties such as 
refraction at a negative angle, an inverse Doppler shift, and a 
backwards oriented Cerenkov radiation cone [1]. Veselago [1] 
also predicted that, under the conditions 
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a parallel slab of this material could focus waves emitted 
from a point source and essentially act as a lens, albeit 
without magnification. These conditions are considered "ideal" 
because they lead to the same impedance and speed of light 
as free space. 

Pendry [4] pointed out that under these conditions incident 
decaying evanescent waves become growing inside the slab. 
The recovery of the evanescent waves allows the image reso- 
lution to be better than the diffraction limit and this material 
could theoretically produce a perfect image of the source, 
which lead Pendry to call the system a "perfect lens." A BW 
material must necessarily be dispersive [1] and therefore all 
this would only work properly at certain frequencies, where the 
conditions of Q are fulfilled. At these frequencies the meta- 
material slab not only compensates for the phase propagation 
of the waves, as a conventional lens does, but also compensates 
for the amplitude decay experienced by evanescent waves. 

This claim of growing evanescent waves and the possibility 
of a perfect lens has aroused much interest and given rise to 
much discussion of its correctness and feasibility, e.g., [5]- 
[12]. Besides that, there is great interest and potential for 
applications in these materials beyond the issue of a perfect 
lens. 

In this paper we study the behavior of an evanescent wave 
interacting with a slab of BW material via simulations. We 
consider the performance and applicability of two numerical 
techniques being used to simulate this interaction. In Sec. |ll| 
we introduce the model system and in Sec. [Ill] the numerical 
schemes. Section IIVI discusses results and Sec. [V] presents 
conclusions. 

II. Model System 

To be able to study the properties of evanescent waves 
interacting with a BW slab, we isolate a single evanescent 
wave by inserting the slab into a parallel plate waveguide and 
use excitation below the cutoff frequency of the waveguide. 
This way the wave in the waveguide is evanescent and its wave 
vector can be chosen through the width of the waveguide and 
the spatial profile of the source. Figure ^ shows a schematic 
of the system. The top and bottom of the waveguide are 
perfectly electrically conducting (PEC) plates. PEC walls also 
exist at the sides of the computational domain but these 
are behind perfectly-matched layers (PML) which essentially 
allow the structure to mimic an open waveguide (but the 
evanescent fields considered here decay rapidly enough that 
the performance of the PML is not a critical component in 
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Fig. 1. Schematic view of the model system. The parallel plate waveguide 
with width w is translationally invariant in the z direction. A sheet current 
source at y = —a and a slab of BW material are located inside the waveguide. 



the simulation). The system is translationally invariant in the 
z direction and can be treated as two-dimensional. A sheet 
current density j at y — —a parallel to the x-z plane acts as 
the source for the fields. As our interest is in the behavior of 
evanescent fields, a single evanescent mode is excited using 
the spatial profile of a sheet current matched to that of the 
lowest order evanescent mode. 



ii^, y,t) = zj{t)S{y + a) sin (j^^'j , 



(2) 



where S is the Dirac delta function and w is the width of the 
waveguide. The source location and the waveguide width are 
shown in Fig.^ and j{t) is the time dependence of the current 
source. (In the steady-state situation we deal with a single or 
few wave vectors. For this reason the concern of Ref. [7] does 
not apply here.) 

Composite materials in which the effective electric permit- 
tivity and the effective magnetic permeability are both negative 
in some frequency interval have experimentally been shown to 
exist [2], [3]. Based on this we assume that both the electric 
permittivity and the magnetic permeability are described by 
Lorentz-type frequency dependences as 
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Here Lo^p and tOmp are the plasma frequency, LOea and are 
the resonance frequency, and 7e and 7^ are the absorption 
parameter of the permittivity and permeability, respectively. 
To simplify the model we choose and fir to coincide, i.e., 
LOep = tomp = LOp, LUeo = (^mO = t^o, and 7e = 7m = 7. In the 
special case of c^o = 7 = 0, the conditions Q are fulfilled at 
the design frequency 



(5) 



We take (|5} to define the relationship between the design 
frequency and the plasma frequency even when luq and 7 are 
nonzero but small. 

Our attention is primarily concerned with the design fre- 
quency since this is the only frequency at which a perfect 
focus could possibly be achieved. At a frequency other than 



the design frequency the rate at which evanescent fields grow 
in the slab is not equal to the rate they decay in free space. 
Thus for a source emitting multiple evanescent fields in front 
of the slab, there is no unique point on the other side of the 
slab at which the fields have obtained the same level as at the 
source (i.e., there is no unique image point). 

We assume that the material has abrupt edges and the 
change in constitutive parameters is instantaneous and collo- 
cated. The case of the change in the constitutive parameters of 
a BW slab not coinciding leads to a shift in the frequency of 
the nonreflecting wave and the bound modes of the slab [13]. 
The fields of the bound modes decay exponentially away from 
the slab. The two bound modes with spatially decaying fields 
inside and outside the slab delimit a frequency interval that 
also contains a nonreflecting wave with exponential spatial 
dependence. 

For the material reaction to the field to be causal, the 
Kramers-Kronig relations [14], [15] show that for any de- 
viation from free-space behavior the imaginary part of the 
permittivity or permeability cannot vanish for all frequencies. 
On the other hand, as far as causality is concerned, it can be 
arbitrarily small. 

III. Numerical Schemes 

The numerical schemes used here are applied directly to 
Maxwell's equations rather than to a wave equation. The first- 
order partial differential equations are discretized in time and 
space and an update algorithm is developed. Specifically the 
behavior of the traditional Yee finite-difference time-domain 
(FDTD) algorithm [16], [17] and the pseudospectral time- 
domain (PSTD) method [18] are studied. In both schemes the 
system is modeled as two dimensional and the open sides 
of the waveguide are modeled by using a causal uniaxial 
perfectly matched layer (PML) absorbing boundary [19]-[21]. 
One significant difference between the two schemes lies in 
the grid location where the fields are sampled. In the PSTD 
scheme the sample points of the fields are collocated while in 
the FDTD scheme they are staggered. 

The plates of the parallel plate waveguide are modeled 
as perfect electric conductors. On the PEC boundary, the 
tangential components of the electric field and the normal 
component of the magnetic field must vanish. The boundary 
conditions of the other components are given by surface charge 
density ps and surface current density J5, 
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Here n is the surface normal unit vector. 

A. Yee Scheme 

In 1966, Yee [16] proposed a discretization scheme for 
Maxwell's equations based on a fully staggered grid, where the 
electric field components are each sampled on one of the edges 
of a cubical unit cell and the magnetic components each on 
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one of the face centers. In this manner the spatial derivatives 
can be approximated as centered differences and one achieves 
second-order accuracy in the spatial step size [17]. 

The PEC boundary conditions are implemented by placing 
the boundary in a plane that contains the tangential electric 
field sampling points and then setting these components to 
zero. In this case, the tangential magnetic field components 
and the normal electric field component are not sampled on 
the boundary and one does not need to enforce any boundary 
condition on them. The normal component of the magnetic 
field is sampled on the boundary and is updated in the 
usual fashion. Its boundary condition is automatically fulfilled 
through its dependence on the tangential electric field. 

Owing to the staggering of the field components, the 
precise location of an interface where both the permittivity 
and permeability change is ambiguous. The case where a 
single material parameter, such as the permittivity, changes has 
been well-studied in the literature, e.g., [22], [23]. It is well 
established that for a tangential field component collocated 
with the interface, the best approach is to use the arithmetic 
average of the material properties to either side. (It can also 
be shown, at least for the case of normal incidence, that if one 
is free to assign the material interface to fall halfway between 
tangential nodes, an abrupt change in the material parameters 
is optimum and, in fact, superior to averaging [24].) For the 
study here, the precise location of the interface is not a primary 
concern. Instead, the important issue is whether the evanescent 
fields exhibit the proper growth within the BW slab (whatever 
its thickness or location). As reported in [25], averaging of 
material parameters has been used for modeling BW slabs with 
the Yee algorithm. (In that work surface waves are present and 
no subwavelength imaging is seen.) 

Figure |2] shows a portion of the computational grid as well 
as the material properties associated with each node — one 
must attribute a permeability to each magnetic field node and 
a permittivity to each electric field node. Figure|2ja) shows the 
case of using abrupt discontinuities in the material parameters 
in the Yee grid while (b) shows the case where averaging is 
used for the tangential node at the interface. Figure|2jc) shows 
the grid for the PSTD simulation which is described next. 
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Fig. 2. Depiction of the way tlie permittivity and permeability cliange at an 
interface. The horizontal bars show the material parameters associated with 
a particular field component (the associated field is indicated by the label 
to the left of the bar). A slice of the grid indicating the node positions is 
given for each of the cases, (a) The Yee grid with "abrupt" boundaiies. The 
interface is assumed to exist between the last Hx which uses permeability 
and the first node which uses permittivity t2- (b) The Yee grid using 
averaging. The interface is assumed to be collocated with the Hx node which 
uses the average permittivity to either side even though the introduction of a 
new material (/^avg) arguably creates an additional interface further to the left, 
(c) The PSTD grid. Because the grid is not staggered, the boundary between 
the two regions is unambiguously assumed to exist between the nodes which 
employ the different material parameters. 



B. Pseudospectml Scheme 

In pseudospectral techniques a discrete Fourier transform 
is used to calculate the spatial derivatives in the discretized 
version of a partial differential equation. Time integration is 
achieved using the same second-order time-stepping approach 
as used in the Yee algorithm. Unlike the Yee algorithm where 
central-differences are used, in PSTD the spatial derivatives 
are calculated at the same location as the field. As shown 
in Fig. Iljc), for Maxwell's equations this means that the 
electric and the magnetic field components are all sampled 
at the same location, i.e., in the center of a cubic unit cell. 
An important problem with this approach is that the use of 
exponential Fourier transforms inherently introduces periodic 
boundary conditions that are often undesired. This problem 
can be avoided by covering the numerical domain boundaries 
with a PML layer such that the waves are absorbed in the 



PML before they can contaminate the simulation through the 
periodic boundaries [18]. This gives one a tool to simulate 
open domains. If the system contains PEC boundaries, the 
problem becomes much more difficult. Nonetheless, in a two- 
dimensional system in transverse electric (TE) polarization, 
i.e., E perpendicular to the plane, which is source-free near the 
boundaries, one can show that the PEC boundary conditions 
(Eqs. |6H9} imply that 

(n- V) (n X H) = 0. (10) 

Here n is the surface normal unit vector and fi • V is the 
derivative in that direction. Conditions and can be 
easily enforced with Fourier sine and cosine transforms if one 
chooses Fourier sine transforms for the tangential electric field 
components and Fourier cosine transforms for the tangential 
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magnetic field components. No spatial derivatives of normal 
components appear and no boundary condition needs to be 
enforced on them. Thus, in a two-dimensional system in TE 
polarization, the PEC boundary conditions are conveniently 
implemented through the choice of Fourier transform. 

The discrete Fourier transform has problems representing a 
delta function, discretized as a Kronecker delta, correctly. For 
this reason it is advantageous to use a spatially smoothed rather 
than a highly localized source [18] to represent the source 
screen of This problem also appears when a sudden change 
in material parameters leads to a sudden change in the fields. In 
the simulations one observes wiggles or ripples on the signal 
that can become significant if the change in field is strong. This 
artifact is the Gibbs phenomenon associated with the Fourier 
transforms. 



C. Dispersion Implementation 

To achieve negative e and negative /i in a material, this 
material must necessarily be dispersive [1]. In our model 
system we chose Lorentz dispersion characteristics ( I3I4> [2], 
[4]. There are several commonly used methods to implement 
frequency dependence in a time-domain algorithm, e.g., the 
auxiliary differential equation (ADE) method [26], the re- 
cursive convolution method [27], [28], and the Z transform 
method [29]-[31]. In this study we use the ADE method, 
the Z transform method, and two frequency approximation 
methods [32]. Here we present the update for the electric 
field; analogous equations hold for the magnetic field. In 
the following, will denote the temporal step size of the 
discretization scheme. 

1) Auxiliary Differential Equation Method: One approach 
to implement a frequency dependent response in a discrete 
time-domain method is to transform the frequency-domain 
constitutive parameters into the time domain and then approx- 
imate the derivatives with differences based on Taylor series 
expansions [26]. The frequency dispersion (|3} at any given 
electric field node inside the material slab is implemented as 



.o^A?+7A, + 2)^-4i^ 



c^^A? - J At + 2] 



4i?"-i - [{ul + ul) A? - 7A, + 2] i?"-^! 



[{Lol + u:l)Al+^At + 2] 



(11) 



where D is the electric flux density which must be stored as 
a separate quantity. The superscripts on the fields denote the 
time step at which the field values are taken. 

2) Z Transform Method: The Z transform is essentially 
a more general version of the discrete Fourier transform. It 
is strictly based on a sampled signal and provides an exact 
transformation between sample domain and Z domain. For fre- 
quency dependencies of the Debye, plasma, and Lorentz type 
and even some types of non-linearities it allows one to derive 
exact update algorithms (i.e., exact in the sampled sense) for 



the sampled signal [29]-[32]. The frequency dispersion of 
is implemented as 

,2 

(12) 



jn ^ 2e-°'^' cos (/3At) /"-^ - g-^^A. jn-2 

+ e-"'^' sin {l3At)E''\ 



(13) 



where a = 7/2, /3 = ^Jojq — 7^/4, and / is an auxiliary 
field. This implementation is applicable for real /3, i.e., when 
the system response function is a damped sine wave. When (3 
is imaginary, the character of the response function changes 
and the Z transform has to be implemented differently [31]. 

3) Frequency Approximation Method: Another common 
method to derive discrete time-domain update equations for the 
frequency dependence of a system is frequency approximation. 
Here one starts with the frequency domain expression for the 
system response and makes the transition to the discretized 
time domain by replacing each occurrence of iu) by — , 
i.e., a linear backward difference where represents a 
shift back in time of one temporal step [31], [32]. When the 
frequency dependence becomes complicated this method can 
make the derivation of update equations much easier than the 
Z transform method. On the other hand, in general, while it 
does not significantly change the computational burden, it in- 
troduces new error A bilinear approximation, iu — > Slirf^' 
gives better accuracy but increases the computational burden 
[31], [32]. The frequency dispersion Q using the linear 
approximation is implemented as 



7At + l) 



n—1 I T71—2 



(7At + 2)/ 



(14) 



/" = {w^A^^;" + (7At + 2) /"-I - 

X [a;2^2^7At + l]"'. (15) 
In the bilinear approximation the electric field is updated using 



2 a2 ipn-l 



E" = {io'^Af + 2^ At + 4) — ~ 2ujlAiE 

+ {2ujlA^t - 8) + KA? - 27 At + 4) /"-^l 

x[{ul + Lol)A^t+'^l^t + ^Y\ (16) 
/" = {wpA^E;" + 2cj2a2£;"-i + cj2a2£;"-2 

- {2u:lAl ~ 8) /"-I - K^A? - 2^ At + 4) /"-^j 

X [wgA2 + 27At+4]"'. (17) 

In (I12t -( ll7t / is used as an auxiliary field. 

IV. Results 

We have performed model calculations for a slab of BW 
material inside a waveguide using the geometry and the 
material described in Sec. |ll| We compare the results for the 
FDTD and the PSTD schemes, using the ADE method for the 
dispersion implementation. We also show a comparison of all 
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four dispersion implementations of Sec. Illl-Cl with the PSTD 
scheme. 

For the FDTD method one commonly uses the Yee scheme 
[16] to derive the update equations for the fields. Because 
this scheme is based on a staggered grid, where all field 
components are sampled at different locations in a unit cell, 
it has inherent difficulties representing a collocated change 
in material parameters. This problem becomes clearly evident 
when modeling an interface between free space and a BW 
material. 

In the PSTD method all the field components are sampled 
at the same location and no transition layer exists. We find 
that, unless the dispersion implementation introduces phase 
error, the bound modes and the nonreflecting wave appear at 
frequencies very close to those found in the continuous world. 

In the simulations the design frequency is 15 GHz and 
the spatial discretization level is 100 points per propagating 
free-space wavelength at the design frequency A/, i.e., Aj. = 
Ay = A//100 = c/(100/rf) = 0.1999mm. Using such a fine 
discretization is typically unnecessary when simulating propa- 
gating waves, where it is more customary to use approximately 
20 points per wavelength, but because evanescent waves have 
rapid amplitude decay and also phase variations (transverse to 
the amplitude decay) that have wavelengths smaller than the 
propagating free-space wavelength, it can be necessary to use 
a finer discretization to model these waves accurately. 

The slab length L is 33Aj,, and the source is located a 
distance a = 3QAy in front of the BW slab. The waveguide 
width w which determines the degree of evanescence was 
chosen to be 32Ax- (For w > \f /2 — 5QAx the waves in the 
waveguide would be propagating.) The material parameters 
were chosen as w„ /(27r) V2fd, t^o/(27r) = 5 MHz, and 
7/ {2tt) — 5 MHz. The values of ujo and 7 are small compared 
to ujp and their influence on the results is weak. Even though 
Q can only be approximately fulfilled, due to the loss and 
the nonzero resonance frequency, we will assume that the 
design frequency is unchanged. The results for the Lorentz 
material are similar to those obtained for an unmagnetized 
plasma, i.e., luq — OHz. The time dependence of the source 
current is a sinusoid with an exponential ramp j{t) = (1 — 
g-t/r-j sm{27rfst), with carrier frequency fs and the ramping 
governed by r. In the simulations we used r = 5000 At. 

The graphs presented in this section show the spatial depen- 
dence of the magnitude of a certain frequency component of 
the temporal Fourier transform of the fields in the y direction. 
The temporal Fourier transform was recorded between 15000 
time steps and when the simulation stopped at a total of 67000 
time steps. The solid vertical lines indicate the extent of the 
BW slab, the dotted vertical line indicates the source location, 
and the arrows show an arbitrarily chosen object location 
and its corresponding image location. Ideally, at the design 
frequency the BW slab compensates exactly for the free-space 
propagation and decay of incident waves corresponding to 
the thickness of the slab. Thus, referring to the geometry of 
Fig. n an object located at y — —b with 6 < L, is imaged at 
y = 2i - 6. 
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Fig. 3. Magnitude of tlie Fourier tansform component at the design 
frequency using tlie abrupt-boundary Yee scheme with the ADE dispersion 
implementation. The lines for the different fields have been normalized and 
offset to allow qualitative comparison. The system parameters are given in 
the text. The solid vertical lines indicate the extent of the BW slab, the dotted 
vertical line indicates the source location and the arrows show an arbitrarily 
chosen object location and its corresponding image location. 



A. Nonreflecting Wave 

In the Yee discretization scheme, all field components are 
sampled at different locations in a unit cell. The permittivity 
and permeability are usually also sampled at these locations. 
A change in material parameters at the interface between two 
materials, therefore, is extended in space and the algorithm 
essentially sees a thin transition layer, where some parameters 
are those of one material and the others are those of the 
other material. As discussed in Sec. IIII-AI the interface may 
be approximated with or without an average of the material 
properties to either side. We identify the approach without 
averaging as the "abrupt" boundary. 

For a BW slab, the effect of such a transition layer between 
the change in e and the change in fi has previously been stud- 
ied using a continuous-space model [13]. Our Yee simulations 
employing abrupt boundaries have transition layers with the 
e of the BW slab and the /i of the surrounding free space, 
similar to the situation in [13]. For the front face of the slab, 
this corresponds to the situation shown in Fig. |2ja) if one 
associates the subscript 1 with free space and the subscript 2 
with the BW material. When averaging is used, the average 
permeability is applied to the nodes assumed to coincide 
with the interfaces. 

To illustrate the effect of the transition layer in the Yee 
scheme, we show in Fig. |3] the fields at the design frequency, 
where the source frequency fs equals fd, using abrupt bound- 
aries and the ADE dispersion implementation. The evanescent 
field grows inside the BW slab but a reflected field is clearly 
present, as can be seen from the change in slope of the fields 
between the source and the y = surface, as well as between 
the two slab surfaces. In Fig. |3] the field amplitudes at the 
image location are significantly different from those at the 
object location. Figure 0] again shows the results at the design 
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Fig. 4. Same as the previous figure except except averaging is used for tlie 
permeability at the interface. 
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Fig. 5. Magnitude of the Fourier transform component at / = 14.83 GHz, 
away from the design frequency = 15 GHz, using the abrupt-boundary Yee 
scheme with the ADE dispersion implementation. The lines for the different 
fields have been normalized and offset to allow qualitative comparison. The 
system parameters are given in the text. 



frequency except now averaging of the permeability is used at 
the interface. As before, the results do not exhibit the desired 
behavior The deep nulls in both these figures occur when 
the total field changes sign because the reflected field is of 
equal magnitude as the field of the source at that point but 
has opposite sign. 

Figure |5] shows the nonreflecting wave in the Yee scheme 
with the abrupt boundary. It occurs at 14.83 GHz, i.e., below 
the design frequency fd of 15 GHz, and was found by scanning 
through the frequency range between the two bound mode 
frequencies. The field is clearly amplified exponentially inside 
the BW slab and decays exponentially outside of it. Upon 
close inspection one finds that the magnitude of the electric 
field at the object location is approximately 28, while it is 
about 26 at the image location. This discrepancy is due to the 



Fig. 6. Magnitude of the Fourier transform component at / = 15.17 GHz 
using the Yee scheme with material averaging and the ADE dispersion 
implementation. The lines for the different fields have been normalized and 
offset to allow qualitative comparison. The system parameters are given in 
the text. 



slower speed of light in the BW material at this frequency. 
Thus, the exponential growth inside the slab is slower than 
the decay outside of it and, though the field is amplified in 
the slab, the amplification does not accurately compensate for 
the decay in free space. The magnitude of this discrepancy as 
well as the frequency at which the nonreflecting wave occurs 
is dependent on the wavevector of the field [13]. Thus, it is 
not possible to define a single corrected image location that 
holds for all wavevectors (this is true of both the propagating 
and evanescent components). 

Figure |6] shows the nonreflecting wave when averaging is 
used for the boundary. Here the nonreflecting wave was found 
to occur at 15.17GHz, i.e., above the design frequency. An 
object of field strength 28 was found to have a corresponding 
field strength at the image point of 33 (i.e., the field was 
overcompensated). 

In [13] it was found that transition layers shift the frequency 
of vanishing reflection coefficient upward from the design 
frequency. In the FDTD simulations we find that for the 
abrupt-boundary Yee scheme the nonreflecting frequency is 
shifted downward from the design frequency. We attribute this 
difference to the effect the discretization has on the transmis- 
sion and reflection coefficient of each interface. (This is in 
addition to introducing transition layers into the system.) When 
averaging is used, the nonreflecting frequency is above the 
design frequency, but the offset between the two frequencies 
is approximately the same as when the abrupt boundary is 
used. Thus there does not appear to be an obvious advantage 
to using averaging when modeling a BW slab since neither 
yields the correct behavior at the design frequency. (However, 
again note that we are not concerned with the specific location 
of the boundary, but rather just the behavior of the evanescent 
fields.) 

In the PSTD scheme all field components are collocated and 
there is no transition layer at a material interface. Using the 
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Fig. 7. Magnitude of the Fourier ti'ansform component at the design 
frequency using the PSTD scheme with the ADE dispersion implementation. 
The lines for the different fields have been normahzed and offset to allow 
qualitative comparison. The system parameters are given in the text. The 
ripple on Hx is a numerical artifact of the Fourier transform due to finite 
spatial resolution. 



ADE method to implement the frequency dispersion of the BW 
material we observe that the nonreflecting wave occurs at the 
design frequency, as shown in Fig.0 Because the wave vector 
magnitudes inside and outside the material are nearly equal 
(with slight discrepancy due to loss and a nonzero resonance 
frequency), the amplification of the evanescent wave inside the 
slab compensates for its free-space decay. The field amplitudes 
at the object location and its image location are within 0.5% 
of each other. 

The examples in Figs. are at a discretization level of 
A//100. Arguably, for either the averaged or abrupt boundary, 
this provides a transition layer thickness of only A//200 in 
the Yee scheme and yet it still produces a significant enough 
frequency shift for the nonreflecting wave that realization of 
a perfect planar lens is essentially impossible. The analysis 
in [13] shows that in the limit of vanishing layer thickness 
the frequency of the nonreflecting wave goes to the design 
frequency in the continuous case. Nevertheless, for all prac- 
tical discretization levels the influence of the transition layer 
remains significant. There have been reports of the Yee scheme 
being used to model evanescent fields in BW materials [33], 
[34], but the level of discretization that had to be used was on 
the order of 566 cells per propagating wavelength [35] which 
may be prohibitive for most applications. 

In general all finite difference schemes suffer numerical dis- 
persion. The PSTD scheme is more computationally expensive 
than the Yee algorithm for the same level of discretization but 
suffers less numerical dispersion. The dispersive properties of 
the PSTD scheme have been discussed in the literature [18] 
(although the dispersion relation given in the literature only 
pertains to the case of an infinite grid and thus only holds 
approximately for finite grids). The dispersive properties for 
evanescent fields in FDTD simulations have also be discussed 
in the literature [36]. The fields being considered here are 



BW Slab 




Position (A^) 

Fig. 8. Magnitude of the Fourier transform component at the design 
frequency using the PSTD scheme. The field is calculated using the ADE, Z 
transform, and the bilinear and the linear frequency approximation methods. 
A reflected wave is clearly visible with the linear frequency approximation. 
The system parameters are given in the text. The lines for the different fields 
have been offset to allow qualitative comparison. 

discretized such that numeric dispersion is not a significant 
concern. Thus the major difference between the schemes 
considered here is the treatment of the interface and one 
cannot attribute the superior results of the PSTD scheme to 
its superior dispersion. 

B. Comparison of Dispersion Implementations 

In Sec. IIII-CI we discussed four different implementations 
of the frequency dependence of the permeability and permit- 
tivity. To compare their performance, we use each one in a 
simulation using the PSTD scheme at the design frequency, 
with the parameters given at the beginning of Sec. IIVI In 
this configuration one expects the reflected wave at the first 
interface {y = 0) to vanish, as shown in Fig. In Fig. |8l 
we show the electric field in these simulations. It is apparent 
that the ADE, the Z transform, and the bilinear frequency 
approximation all give good results. No reflected wave is 
visible and upon inspection one finds that the magnitudes at 
the object and image location are essentially equal. When 
plotted without offset the three curves fall on top of one 
another The linear frequency approximation on the other hand 
does show a significant amplitude for the reflected wave. The 
field at the image location is also much weaker than at the 
object location. This effect is probably due to the numerical 
phase error inherent in this method [32]. This shows that the 
use of the linear frequency approximation is problematic in 
simulations involving BW materials, while many of the other 
common dispersive material implementations work well. 

V. Conclusion 

We have studied the finite-difference time-domain and the 
pseudospectral time-domain methods in the context of evanes- 
cent waves in the presence of a backward-wave material. Our 
simulations verify that the predicted growth of the evanescent 



8 



ACCEPTED BY ffiEE TRANSACTIONS ON ANTENNAS AND PROPAGATION FOR JANUARY 2005 



field inside the BW material, which is vital for imaging beyond 
the diffraction limit, does occur. We found that the Yee FDTD 
method suffers from a transition layer effect due to the inherent 
staggered grid, even when using an averaging technique at 
the material interface, while the collocated grid used in the 
PSTD method avoids this numerical artifact. The choice of 
numerical technique to implement the frequency dependence 
of the permittivity and the permeabiUty is also important for 
the correct modeling of the BW material. 

Thus, we have shown that the PSTD method in concert with 
the ADE technique, the Z transform technique, or the biUnear 
frequency approximation, allows one to accurately model the 
interaction of a BW slab with an evanescent wave. 
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